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Abstract. We study the gas-surface dynamics of O2 at Ag(lll) with the 
particular objective to unravel whether electronic non-adiabatic effects are 
contributing to the experimentally established inertness of the surface with respect 
to oxygen uptake. We employ a first-principles divide and conquer approach based 
on an extensive density-functional theory mapping of the adiabatic potential 
energy surface (PES) along the six O2 molecular degrees of freedom. Neural 
networks are subsequently used to interpolate this grid data to a continuous 
representation. The low computational cost with which forces are available from 
this PES representation allows then for a sufficiently large number of molecular 
dynamics trajectories to quantitatively determine the very low initial dissociative 
sticking coefficient at this surface. Already these adiabatic calculations yield 
dissociation probabilities close to the scattered experimental data. Our analysis 
shows that this low reactivity is governed by large energy barriers in excess 
of 1.1 eV very close to the surface. Unfortunately, these adiabatic PES 
characteristics render the dissociative sticking a rather insensitive quantity 
with respect to a potential spin or charge non-adiabaticity in the 02-Ag(lll) 
interaction. We correspondingly attribute the remaining deviations between 
the computed and measured dissociation probabilities primarily to unresolved 
experimental issues with respect to surface imperfections. 
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1. Introduction 

A predictive materials science modelling based on microscopic understanding requires 
a thorough knowledge of all underlying elementary processes at the atomic scale. 
The (dissociative) adsorption of individual oxygen molecules at metal surfaces is such 
an elementary process that is of crucial relevance for a wide range of applications, 
with heterogeneous catalysis forming just one prominent example. Recent work on 
the dissociative adsorption of O2 at Al(lll) has severely challenged the prevalent 
understanding of this process [TJ [U |3J [3]. Significant spin non-adiabatic effects, 
i.e. a coupled electronic-nuclear motion beyond the standard Born-Oppenheimer 
approximation, have been proposed as rationalization for the experimentally observed 
low initial adsorption probability. This behaviour has been associated to the low 
density-of-states (DOS) at the Fermi-level of the simple metal surface, which prevents 
an efficient spin quenching through the tunnelling of electrons between substrate and 
adsorbate. This hinders the transition from the spin-triplet ground state of gas-phase 
O2 into the singlet state upon adsorption at the metal surface. The low sticking 
coefficient can then be seen as a result of similar spin-selection rules that are well 
known from gas-phase chemistry. 

In particular from the point of view of heterogeneous catalysis, interest naturally 
shifts to low-index noble metal surfaces. Only recently, a spin- asymmetry in electron- 
hole pair excitation spectra for O2 at Pd(100) have further fortified the previously 
proposed DOS-related mechanism [5]. On the other hand, in contrast to palladium, 
the DOS at the Fermi-level is equally low for aluminium and coinage metals. This 
raises the question whether similar spin non-adiabatic effects significantly contribute to 
their established low reactivity with respect to O2 dissociation [6ll7]l8]li^lT0 l lTT | [ll ^ [T3 ] . 
Here, we aim to contribute to this with a first-principles based modelling of the O2 
gas-surface dynamics specifically at the Ag(lll) surface, for which particularly low 
adsorption probabilities have been reported experimentally [101 HI] • Earlier simplified 
models based on low dimensional analytic potentials have already tried to understand 
the nature of such low reactivity in terms of charge transfer limitations from the 
Ag(lll) surface to the impinging O2 molecule [HI [16]. Notwithstanding, on Ag(100) 
recent molecular dynamics (MD) calculations based on an ah initio six-dimensional 
potential energy surface (PES) showed that the comparable absence of O2 dissociation 
on that surface can be fully explained within an adiabatic picture |17j . There, the 
system is simply characterized by large dissociation barriers of about 1.1 eV that 
appear when the molecule is close to the surface. This motivates a similar ah initio 
investigation also for Ag(lll), with the objective to settle the question up to what 
extent the low reactivity on the Ag(lll) surface is really a signature of electronic 
non-adiabaticity or only the result of large dissociation barriers as found for Ag(100). 

The overall structure of the paper is as follows. The next section provides a 
detailed description of the different steps employed to model the dynamics of O2 on 
Ag(lll) within a divide-and-conquer approach. Section [3] commences with a brief 
review of the existing experimental data, emphasizing controversies with respect to 
the absolute order of magnitude of dissociative sticking for incidence energies Ei 
below 1 eV. We then proceed with a discussion of the quality and relevant static 
properties of the computed adiabatic PES, revealing as a central characteristic barriers 
to dissociation larger than 1.1 eV at close distances to the surface. The ensuing 
classical trajectory calculations demonstrate that these barriers lead to dissociation 
probabilities that already fall within the large error bars set by the scattered 
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experimental data. With heights that are larger than the gas-phase O2 spin singlet- 
triplet gap (0.98 eV [IS]) and the 2 electron affinity (0.44 eV [TSJ) these barriers 
furthermore completely mask any signs of a potential spin or charge non-adiabaticity 
in the 02-Ag(lll) interaction in the dissociative sticking coefficient. The remaining 
deviations of the computed sticking curve to the measurements for Ei < 1 eV are thus 
more likely the result of the unresolved issues in the experimental data, e.g. with 
respect to surface imperfections as discussed by the experimentalists [H] [19] . 

2. Methodology 

Methodologically, the very low dissociation that characterizes this system makes the 
use of on-the-fly ah initio molecular dynamics (AIMD) to determine the sticking 
coefficient |20) computationally unaffordable. Instead, our simulations are based on a 
divide-and-conquer approach [2T1E21 23J, in which an accurate PES is first constructed 
from first-principles, and next, classical MD calculations are performed using a 
continuous representation of this PES to describe the molecule-surface interaction. 
Once the former is obtained, the latter come at negligible computational cost. In 
contrast to AIMD, the divide-and-conquer approach thus enables the calculation of 
a large number of trajectories. Even for very small sticking probabilities and various 
initial conditions of molecules in the impinging molecular beam (Ei, incidence angle 
etc.), sufficient statistical sampling can thus be performed. Additionally, detailed 
information on the PES properties, which allows for a better rationalization of the 
factors ruling the gas-surface dynamics, provide another advantage. Simulations based 
on this methodology have been able to reproduce not only reactive probabilities, 
but also more subtle quantities like the rovibrational state distributions of scattered 
molecules [231125]. 

The key ingredient of this impressive success that ensures reliable quantitative 
results is the use of multi-dimensional ah initio-b&sed PESs. In this respect, several 
interpolation methods to obtain continuous representations in the context of gas- 
surface dynamics have been developed in the last years, reaching energy precisions on 
the order of meV in parts most relevant for the former [26, 27, 28][2JJ1[3D]. In the present 
work, the PES interpolation is based on neural networks (NN) [28] [29] [30] that we 
adapt to the specific symmetry of the 02/Ag(lll) system as explained in section I2T21 
In order to keep the multi-dimensionality of the 02/Ag(lll) PES in a tractable 
limit, we follow common practice and adopt the frozen surface approximation. This 
approximation has been successfully applied to simulate the reactivity of various 
diatomic molecules (H 2 , N 2 , O2, . . . ) on metal surfaces Q] [24j [31] [32] [33]. One 
rationalization for this approximation has been that the critical bond activation, that 
will irreversibly lead to dissociation or the reflection back into the gas-phase, often 
occurs at short interaction times, for which substrate mobility is not yet a critical issue. 
Obviously, the actual dynamics that proceeds after such a critical activation and closer 
to the surface is then described wrongly - the more, the heavier the adsorbate. But this 
has less consequences for the sticking probability as long as only the distinction "direct 
dissociative adsorption" or "immediate back-scattering" matters. The situation is less 
clear in case of molecular trapping that typically involves long interaction times. We 
return to this problem in section 12.31 below. 
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Figure 1. (Colour online) Coordinate system used in the mapping of the 6D 
Oa/Agflll) PES. Surface Ag atoms are shown as large grey spheres, and the O2 
molecule as small red spheres. The origin of the coordinate system is placed on a 
Ag surface atom. 

2.1. Adiabatic PES from DFT calculations 

Within the frozen-surface approximation the remaining six-dimensional PES only 
depends on the degrees of freedom of the O2 molecule, given by the cartesian 
coordinates Ra = (Xa, Ya,Za) and Rb — Yb, Zb) of the two constituent oxygen 
atoms, A and B, respectively. In the following, another very common equivalent 
coordinate system as depicted in figure [T] is made use of: Molecular configurations are 
conveniently described by the centre of mass of the molecule (X, Y, Z) , the interatomic 
distance d, and the orientation relative to the surface, defined by the polar and azimuth 
angles $ and tp, respectively. The origin is placed on top of a silver atom within the 
surface plane, and ip) — (90°, 0°) corresponds to the molecular axis oriented parallel 
to the [110] direction. To construct the 6D PES, we interpolate with the NN a large 
set of energies that have been first calculated with density-functional theory (DFT) for 
different positions and orientations of the molecule over the surface. Note that the role 
of subsurface oxygen [MJ [35] is not considered in this work because we are interested 
in the measured initial sticking coefficient, i.e. in the dissociation probabilities that 
are measured at very low coverages (about 0.1 ML or less). We have performed more 
than 7000 spin-polarized DFT total energy calculations, starting with so-called elbow 
plots, i.e. (Z, d) cuts through the 6D PES, at high symmetry sites of the surface. 
Additional irregularly spaced molecular configurations, typically in form of "scans" 
both in lateral (X, Y) and angular (1?, tp) dimensions, were added iteratively as needed 
during the neural network interpolation [vide infra). 

An energy cut-off of 400 eV in the plane- wave basis set together with ultra-soft 
pseudopotentials as contained in the default old Cambridge library of the CASTEP 
code [35] have been used to describe the interaction between electrons and nuclei, 
relying on the exchange-correlation functional due to Perdew, Burke and Ernzerhof 
(PBE) [37]. Based on these computational settings, the computed bulk Ag lattice 
constant and bulk modulus compare well to experimental values [35] ( a^p BE = 4.14A 
and B^| BE = 97.1 GPa vs. a^ xp = 4.07A and s£f xp = 109 GPa, respectively), and 



Non-adiabatic effects during the dissociative adsorption of O2 at Ag(lll)? 



5 



also the free 2 bond-length is with d° p BE = 1.24 A well described (d | = 1.21,A 
[T5]V Notwithstanding, similar to other frozen-core calculations, the 2 gas- phase 
binding energy is underestimated compared to accurate full-potential PBE results [39 . 
With 5.64 eV this fortuitously yields a value close to the experimental gas- phase 
data 00] (5.12 eV). For the binding energies addressed below this is not to be 
of concern, as this underestimation largely cancels in the total energy differences 
constituting the latter, and a potentially introduced constant shift does not affect 
the dynamics on the interpolated PES. 

The Ag(lll) surface is modelled with a periodic five-layer slab separated by 12 A 
of vacuum. In order to sufficiently suppress spurious interactions of the impinging 
O2 molecule with its periodic images in the supercell geometry we are forced to 
use a (3 x 3) surface unit cell. The irreducible wedge of the first Brillouin zone 
corresponding to this large unit cell is then sampled with a (4 x 4 x 1) Monkhorst- 
Pack grid of special fc-points. Systematic convergence tests focusing on the adsorption 
of atomic and molecular oxygen at high-symmetry sites confirm that the binding 
energies are numerically converged to within 50 meV at these settings. For the clean 
Ag(lll) surface we reproduce the established marginal 1% outward relaxation of the 
first interlayer distance [41 . Within the frozen-surface approximation we neglect this 
and compute the 02/Ag(lll) PES above a bulk-truncated slab. As discussed in detail 
by Behler et al. for O2 on Al(lll) [TJ [2], also in this adiabatic PES the spin of the 
O2 molecule is gradually quenched from the triplet ground-state of O2 at large Z 
distances to the spin unpolarized state close to the surface. We use the configuration 
where the O2 molecule is located in the middle of the vacuum region as zero reference. 
This corresponds to the full decoupling limit between O2 and Ag(lll), since at these 
distances, unlike at Al(lll) [UH], there is no spurious charge transfer to the molecule. 
In the employed sign convention negative PES values indicate exothermicity with 
respect to this limit, while positive values indicate endothermicity. 

2.2. Interpolation by symmetry-adapted neural networks 

Neural networks represent a very flexible class of composed functions that in principle 
allows to approximate (physical) potential energy surfaces to arbitrary accuracy 
[421 143] . However, no physical properties of the latter are built therein a priori. 
The use of NNs in the context of gas-surface dynamics, pioneered by Lorenz et 
al. [28], has been described in detail before [29j [30] . Therefore, we only give a 
brief account of the novel aspects of how symmetry is taken into account for the 
present system. Instead of using the "straightforward physical" coordinates depicted 
in figure [T] the neural network interpolation is performed on a properly transformed 
set, which has been termed symmetry functions before |30) . This way, two molecular 
configurations that are equivalent by symmetry are enforced to have same energy 
on the NN-interpolated PES by construction. By the rigid incorporation of these 
important physical properties of a specific system, the interpolation problem is thus 
simplified considerably (i.e. better quality with less input data). The construction of 
such a coordinate transformation is a rather complicated task, as it has to capture 
both periodic and point group symmetry on high-symmetry sites of the surface 
correctly. This implies certain demands for the behaviour under a certain (discrete 
but infinite) set of transformations including the molecular centre of mass position on 
the surface (X, Y) and the molecular orientation ($, ip), but not the distance from the 
surface Z and the internuclear distance d. Artificial symmetries have been introduced 
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Figure 2. (Colour online) Structure and sixfold symmetry of the "(111)" surface 
employed in the symmetry-adapted NN approach (see text). The primitive lattice 
vectors ai and ai are shown as arrows, high symmetry sites are marked by 
different symbols, and the irreducible wedge is indicated by the dark gray triangle. 
Distances are conveniently described in units of the surface lattice constant a. 



unintentionally in past attempts and only recently been found to have a significant 
effect on the calculated sticking probability [H]. 

For the present system, symmetries resulting from the combination of the 
threefold symmetry of the (111) surface of the closed-packed fee metal silver and 
the exchange symmetry of the homonuclear O2 molecule can be further increased by 
exploiting the near degeneracy of the O2 interaction with the fee and hep hollow sites: 
For equivalent molecular configurations at these two sites, we observe a statistical 
relation in our DFT data set with a root mean square error (RMSE) of less than 
5 meV and a mean absolute deviation of less than 3 meV. Similar results have been 
obtained before for H2 on Pt(lll) [45]. For the present work we will neglect this 
small difference and thus treat the fee and hep sites as equivalent. This increases the 
symmetry of the surface from three- to sixfold, as illustrated by the high-symmetry 
points and irreducible wedge in figure [5J together with the primitive lattice vectors 
ai and &2- Inspired by the symmetry functions obtained for a fcc(lll) surface with 
threefold symmetry [30], the following set of symmetry-adapted coordinates Qi have 
been designed and carefully verified not to yield any artificial symmetries as "side 
effects" . A detailed derivation and their simple generalization to other common low- 
index surfaces of both close-packed fee and hep metals will be detailed elsewhere [IB] : 
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3! and g 2 are functions R 2 — >• K with the periodicity given by the surface lattice 
vectors ai and a 2 . They have been constructed such that points r = (x,y) in the 
surface plane which are equivalent by the sixfold symmetry depicted in figure [5] are 
mapped onto the same pair of function values (gi(x,y), g 2 (x,y)) [46] : 
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a is the surface lattice constant given by l/V2a^p BE for the present system. The 
building blocks of g\ and g 2 are low order terms of the corresponding Fourier series as 
indicated by the respective reciprocal lattice vectors 



G nm = nbi + mb 2 . 



(2.3) 



Here b\ and b 2 are the primitive vectors of the reciprocal lattice given by the defining 
relation 



ai ■ bj = 2wSij i,j £ {1, 2} 



(2.4) 



3 ^ 

and the crystallographic convention to indicate a negative direction by a bar over the 
corresponding number is adopted. 

Along the lines of previous work |30j . neural network training has been 
simplified by limiting highly repulsive energies to < 5 eV, which is still much larger 
than the energy range of interest for the ensuing molecular dynamics trajectories 
discussed below. Using the adaptive extended Kalman filter algorithm including the 
modifications described in (35] for training, the sensitivity of the NN is further focused 
on the most relevant parts of the PES by assigning training weights 
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Figure 3. (Colour online) Illustration of the interpolation quality obtained with 
the best NN fit for the present system, (a) Representative (d, Z) (elbow) plot 
corresponding to a molecule with X = 0.5 a, Y = \/3/6a and ■& = 0°, where a 
is the surface lattice constant. The size of the black circles shows the maximum 
error of the NN interpolation compared to the DFT input at the corresponding 
configurations, (b) Comparison of the NN interpolated values with all available 
corresponding DFT energies. Chosen to match the estimated accuracy of the 
latter (see text), the grey lines at ±50meV serve as a guide to the eye for the 
interpolation error. The resulting maximum absolute deviation (MAD) and root 
mean square error (RMSE) are given. 



to the input data. More than forty different NN topologies have been tried. In each 
case, the interpolation quality has been carefully monitored by plotting numerous two- 
dimensional (d, Z), (X, Y) and (•&, ip) cuts through the such obtained PESs, including 
information about the interpolation errors of individual data points. An example for 
a corresponding elbow is shown in figure [3] (a). With increasing knowledge gained 
about the PES and dynamics thereon, we decided for a {9-16-16-16-1 tttl} topology 
of the NN. In this commonly used notation [2, 28, 23130], the first number refers to 
the nine symmetry-adapted coordinates (Qi, ■ ■ ■ , Q9) in the input layer. The three 
numbers in the middle denote the number of nodes in the hidden layers, and the last 
number indicates that there is a single node yielding the PES value in the output 
layer. Like in previous work [2] I29j . a hyperbolic tangent and a linear function 
serve as activation functions in the hidden and output layer and are denoted by t 
and I, respectively. Figure 131(b) summarizes the high quality of the finally obtained 
interpolation. For statistically significant subsets of the data points in the dynamically 
relevant entrance channel and chemisorption well (vide infra), the RMSE is smaller 
than 26meV. Considering points from the training and test sets separately, RMSEs 
for the latter are larger by only a few meV - thus indicating good predictive properties 
for arbitrary molecular configurations to be encountered during the dynamics. 

2.3. Molecular dynamics simulations 

On the thus obtained continuous NN-representation of the adiabatic PES we perform 
trajectory calculations using a Bulirsh-Stoer integrator [47] for the classical equations 
of motion in Cartesian coordinates. Trajectories start with the O2 molecule at its 
calculated gas-phase bond length and at Z = 7 A above the surface, neglecting its 
initial zero point energy. A conventional Monte Carlo procedure is used to sample all 



Non-adiabatic effects during the dissociative adsorption of O2 at Ag(lll)? 9 

possible initial O2 orientations (1?, ip) and lateral positions (X,Y). For each incident 
energy the sticking coefficient is obtained based on 10 7 different trajectories at 
perpendicular incidence^ Such a large number of trajectories can only be integrated 
because the cost to compute forces from the NN PES-representation is negligible 
compared to the DFT calculations. As outcome of each individual trajectory we 
consider the following possibilities: 

(i) The trajectory is assumed to yield a dissociation event when the O2 bond length 
d is stretched beyond 2.4 A, i.e. twice the gas-phase distance, and the associated 
velocity d is still positive. This is a rather conservative criterion and we have 
verified that some variation of the critical bond distance has no effect on the 
obtained sticking curves. 

(ii) The molecule is assumed to be reflected when its centre of mass reaches the initial 
starting distance of 7 A with a positive Z- velocity. The results were equally found 
to be insensitive to reasonable variations of this criterion. 

(iii) Molecular trapping, when the molecule is neither dissociated nor reflected after 
24 ps. A detailed analysis of the obtained trajectories shows that this time span 
is well separated from those of reflection or dissociation events, which happen 
in the order of 2-3 ps, such that the precise value of 24 ps has little influence 
on the number of determined trapping events. Although one may expect that 
the longer the molecule spends close to the surface, the more probable would 
be energy dissipation into either electronic or phononic excitations that in turn 
would prevent these trapped molecules to eventually escape from the surface, no 
final statement can be done in such cases within the frozen surface approximation 
(vide infra). 

3. Results and discussion 

3.1. Insight from experiment 

With silver as the almost unique industrial catalyst employed in the selective 
epoxidation of ethylene (48) 09], the interaction of O2 with silver surfaces has already 
been extensively investigated during the last decades in an attempt to understand 
the fundamentals behind such exceptional catalytic functionality (e.g. see [50] and 
references therein). The richness and complexity of the O2 interaction with silver is 
reflected in a variety of adsorbed species that have been identified even at just the 
flat low-index surfaces Ag(100), Ag(110), and Ag(lll), when varying the exposure 
conditions such as surface temperature, oxygen coverage, gas temperature and gas 
pressure (48J EU [52J [53l [54] . In ultra- high vacuum physisorption states are observed 
at lowest temperatures, while molecularly chemisorbed states prevail up to T s ~ 150 K. 
Dissociative adsorption appears for higher crystal temperatures; in this case subsurface 
oxygen can be also important for coverages where silver oxides are about to be formed 
[34] . All these studies showed that the abundance of each state also depends on the 
surface structure, with the close-packed Ag(lll) surface being the most inert towards 
thermal oxygen adsorption [5T] . 

All molecular beam experiments performed on the three low-index faces show 
that, if present, both the molecular and the dissociative initial adsorption probabilities 

\ Results of calculations for various non-perpendicular incidence angles focusing on scattering will 
be reported in a forthcoming publication. 
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are always activated [51 [5J M HO] • Despite this consensus, there remain unresolved 
experimental discrepancies on the Ag(lll) surface regarding the actual order of 
magnitude of the dissociation probability at normal incidence. Using a molecular 
beam, Raukema et al. reported initial dissociation probabilities that increase from 
10~ 5 to 10 -3 as the incidence energy Ei increases in the range 0.8-1.7 eV [10]. In 
contrast, similar experiments performed with a molecular beam of 0.098 < E% < 
0.8 eV, estimated that the dissociation probability on the clean Ag(lll) surface should 
be lower than 9 x 10~ 7 at least in that Ei-range [14] . This has been related [14] to the 
particular propensity of the Ag(lll) surface towards surface imperfections, notably 
steps, or towards residual gas contamination at lower surface temperatures. 

In this situation the main purpose of our modelling is evidently not so much to 
reach full quantitative agreement with experiment, but rather to establish the order 
of magnitude for the initial dissociation probability within an adiabatic picture. As 
we will show below this is already sufficient to reach important conclusions concerning 
the role of possible non-adiabatic effects. In addition it contributes to the discussion 
of possible sources for the differences between the experimental data sets of |T0l [14] 
and will provide a valuable basis for future measurements that aim to overcome the 
present scatter in the experimental data. 

3.2. Adiabatic PES properties 

Overall, the obtained 02/Ag(lll) PES is rather repulsive. A systematic analysis of 
(d, Z) cuts for different molecular configurations shows that at intermediate distances 
from the surface {Z 2 A) the potential energy is already above one eV in almost 
all cases. In this respect, the (d, Z) cuts of figure [3] (a) and figure [4] (a) are rather 
representative. As expected, the repulsion is strongest over the top sites. Already at 
distances Z > 3 A far from the surface the emerging corrugation gives rise to a very 
shallow well with a depth of the order of 20-30 meV. As the employed semi-local PBE 
functional does not account for van der Waals interactions, this is unlikely a proper 
representation of the experimentally reported physisorption state, but instead rather 
an artefact of the NN interpolation. As we will discuss below, this shallow (and likely 
spurious) well only plays a role for the 02/Ag(lll) gas-surface dynamics at lowest 
incidence energies and its effects on the dissociation probability are easily separated. 

Much more relevant for the dynamics is instead a second well that we determine 
at Z = 2.3 A, and for a molecular configuration where the O2 centre of mass is 
situated above a bridge site and with the molecular axis pointing to the two closest 
Ag atoms. Figure 2] (a) shows the elbow plot corresponding to this configuration. The 
computed internuclear distance is with d — 1.28 A only slightly stretched compared 
to the gas-phase value (d — 1.24 A). This indicates a rather weak interaction that we 
find equally reflected in the computed low binding energy of — 40meV and the modest 
shift of the O2 stretch vibration from the computed 192 meV gas-phase value down to 
145 meV. These properties are only little affected by the frozen surface approximation 
and the NN interpolation. After a direct DFT calculation of this state that also 
allows for surface relaxation the most notable change is the increase of the binding 
energy to — 70meV. The latter value is somewhat lower than the — 0.17eV binding 
energy determined before by Xu et al. for the same O2 molecular adsorption state 
[35] . Considering the highly similar computational setup employed in the latter study 
we attribute this 0.1 eV difference primarily to their use of a smaller (2 x 2) surface 
unit-cell. 
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d (A) d (A) 

Figure 4. (Colour online) Elbow plots of two relevant configurations of the 
C>2/Ag(lll) PES: (a) The chemisorption state in which the molecule is located 
over bridge and with (1? = 90°; ip = 0°). (b) Configuration for the minimum 
energy barrier to dissociation in which the molecule is located over bridge with 
(tf = 90°; <p = 90°). 



It is tempting to identify this molecular state with the chemisorption state 
reported in the experimental literature. For this, one has to recognize though that 
(regardless of the small scatter) both theoretical binding energies are at variance 
with the experimental estimate for the chemisorption state of ~ 0.4 eV from thermal 
desorption spectroscopy (TDS) [ST]. This difference is noteworthy as preceding 
equivalent DFT calculations performed for 2 at Ag(100) [17] and Ag(110) [55j [56] 
reproduce the approximately — 0.4eV chemisorption energies measured at these facets 
[48] [57] , We also confirmed that the here employed computational setup fully 
reproduces the O2 binding energy reported by Alducin et al. at Ag(100) |17j . 
This suggests that both the employed PW91 and PBE functionals, which are very 
closely related, are capable of describing the 02-Ag interaction rather well. Similarly 
problematic is the observation that none of the experimentally reported vibrational 
frequencies [SSJ [SJ5] comes close to the 145 meV that we compute for the molecular 
state. Instead much lower frequencies below 100 meV are measured, which, however, 
have already been assigned to hydroxyl, water groups or imperfections |58l I59j . 
Furthermore, as shown by Xu et al. [35. also sub-surface oxygen can strongly stabilize 
chemisorbed O2 at Ag(lll). In this situation we believe that the low computed binding 
energy of the molecular state rather supports the interpretation that the chemisorption 
state prepared in the experiments by Raukema, Campbell and others [101 119[ I51| 153] 
is not a property of the pristine Ag(lll) surface, but reflects the propensity of this 
surface towards defects or residual gas contamination, and we will return to this point 
below when discussing the dissociative sticking probability. Within the theory-theory 
comparison of O2 at the three low-index surfaces, the lower binding energy computed 
here at Ag(lll) as compared to Ag(100) [T7] and Ag(110) [5S][Sj)] is instead perfectly 
consistent with the expected chemical inertness of the close-packed surface. 

As a final important characteristic of the adiabatic PES we also determined the 
minimum energy path (MEP) that leads out of the molecular well to dissociation. 
This path goes along the (d, Z) PES cut shown in figure [4] (b), in which the molecule 
is located over a bridge site with its axis oriented along the direction joining the closest 
fcc/hcp sites. The transition state corresponds to an activation barrier of 1.1 eV, which 
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Figure 5. (Colour online) Initial sticking coefficient So(Ei) of O2 on Ag(lll) as 
a function of the incident energy Ei under normal incidence conditions. Shown 
are the computed probabilities for direct dissociation (black full circles), and for 
molecular trapping in the low (red open circles) and high (green open circles) 
energy range (see text for details). The blue dash-dotted line is included to 
illustrate the total dissociation under the assumption that all of the latter events 
eventually lead to dissociation (i.e. sum of the green open circles and the black 
full circles). The black dashed line illustrates the effect of the finite energetic 
width of the experimental beam on the simulated direct dissociative sticking curve 
according to l|3.1|l . The triangle data points reproduce the experimental data of 
ITOl obtained at T s = 400 K. 



is in very good agreement to the value reported before by Xu et al. [35] . It is located 
at2« 1.6 A, i.e. at a distance closer to the surface {Z ~ 2 — 3 A) than where one 
would expect barriers due to spin or charge non-adiabaticity (vide infra). 



3.3. Dissociative sticking coefficient 

We compile the results of our adiabatic sticking coefficient simulations for normal 
incidence in figure [5] Two different regions can clearly be discerned. At 
incidence energies below 0.4 eV there is an appreciable amount of molecular trapping 
events. This component rapidly decreases with increasing Ei, as characteristic for 
physisorption. Note that the non-monotonous scatter of this curve around Ei ~ 
0.2— 0.4 eV is already no longer statistically relevant. For the employed 10 7 trajectories 
the corresponding probabilities around 10 -7 mean that we have observed either one or 
two molecular trapping events. Analysis of the trajectories in this entire low-energy 
trapping branch reveals that these events are all due to dynamical trapping in the 
shallow PES well located at distances Z > 3 A. As we are rather sceptical about the 
meaning of this shallow well we would not put too much emphasis on this low-energy 
results, even though the order of magnitude of this trapping and the incidence energy 
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range fit rather well to those of the physisorption trapping discussed by Raukema and 
Kleyn [501 HH] ■ We note however that a proper description of any kind of molecular 
adsorption requires to go beyond the frozen surface approximation by including surface 
temperature effects and energy dissipation channels in the simulations. In any case, 
this lowest-energy trapping part is well separated from the sticking properties at 
incidence energies higher than 0.4 eV, and for the latter part the shallow PES well 
plays no further role. 

This other part that we will exclusively discuss henceforth starts at incidence 
energies above 1.1 eV and is characterized by a steeply increasing dissociative sticking 
probability. In contrast, in the intermediate range 0.4 eV< Ei <l.leV we have not 
observed a single dissociation event in the total of 10 7 trajectories performed for every 
incidence energy. From our calculations we can thus put an upper bound to the real 
dissociative sticking probability at these incidence energies of 10 -7 . Above 1.1 eV this 
probability increases steeply and seems to saturate above 10 -2 at the highest incidence 
energies considered. Analysis of the trajectories shows that the dissociation follows a 
rather direct mechanism ruled by the large energy barriers existing at distances below 
2 A from the surface. Almost all of the dissociating molecules follow paths very close 
to the MEP shown in figure BJb) above. This explains the coincidence of the onset of 
sticking above 1.1 eV with the activation energy along this path and rationalizes the 
overall very low sticking probability at this surface in terms of the narrow dissociation 
channel and thus reduced configurational space. 

If we compare these sticking results with the experimental data from Raukema et 
al. [10] that is also included in figure [5] we first notice a gratifying overall agreement 
in the sense that the computed adsorption probabilities are in the right (low) order 
of magnitude and above threshold show the correct monotonous increase up to the 
highest incidence energies considered. However, quantitatively, there are notable 
deviations. This concerns foremost a much steeper slope of the theoretical curve 
especially for Ei < 1.3 eV. As a consequence, there is still appreciable sticking of the 
order 10" 5 - 10" 4 at 0.8-0.9 eV incidence energies in the experiments, whereas the 
theoretical curve drops steeply to below 10~ 7 already at Ei — 1.1 eV. 

One reason that immediately comes to mind to generally rationalize such 
discrepancies is the approximate nature of the employed xc functional. In the present 
case, this is unlikely to explain all of the deviations though. As described before, 
almost all of the direct dissociation goes along paths close to the MEP shown in 
figure 0Jb). If the xc functional messed up the MEP barrier height, this would thus - 
to zeroth order - predominantly have resulted in a rigidly shifted theoretical sticking 
curve to too high or too low Ei, and less to a wrong slope. If the employed PBE 
functional e.g. overestimated the true MEP barrier height, a "better" functional 
would shift the curve to lower Ei. This would then lead to a better agreement at 
the onset region. However, simultaneously it would increase the difference to the 
experimental curve at the highest incidence energies, with the theoretical curve then 
seriously overshooting the experimental data there. 

From this reasoning we thus rather suspect other factors to stand behind the 
disagreement between the theoretical and experimental sticking curve. The influence 
of vibrationally excited O2 in the beam can be neglected, according to the low 
populations (less than 10 -3 for v = 3 and than 10~ 2 for v = 2) that have been 
estimated by Raukema et al. based on the nozzle temperature [TO]. Another factor 
that could be important in particular in light of the abrupt initial increase of the 
theoretical sticking curve is the finite energetic width of the experimental beam [62] ■ 
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To mimic the effect of the latter in the theoretical data we recompute the curve simply 
as a convolution with an energy distribution function /s,.e (-E) for the beam: 



This equation is evaluated numerically, and Sq(E') describes the corresponding curve 
plotted in figure [5] (i.e. a piecewise linear interpolation based on the explicitly 
calculated points). fp : E {E) is taken to be a shifted Maxwell-Boltzmann distribution 



Formulated in the velocity domain, it is known to properly describe supersonic 
molecular beams [53] and has also been used by Raukema for the evaluation of time-of- 
flight (TOF) experiments for the present system [53]. N(ft,E ) = J™ dE' fp, Eo (E') 
is a normalization factor, and the width parameter ft is chosen to reproduce the 
experimental full width at half maximum (FWHM) of 0.2 eV in the energy range of 
interest [60] . As seen in figure [5] this correction leads to a smoother theoretical curve 
that has its onset with 0.9 eV now at slightly lower incidence energies. Nevertheless, 
it does not help to reconcile experiment and theory especially at these Ei close to 
threshold, where the experiment by Raukema et al. [10] yields a dissociative sticking 
that is more than two orders of magnitude higher than the computed one. Our results 
are instead more compatible to those of [2], which estimated a So(Ei) below 9 x 10~ 7 
for Ei = 0.8 eV. 

Another factor that should be kept in mind is the finite surface temperature T s 
in the experiments that is not accounted for at all in the present theoretical approach. 
For the normal incidence data shown in figure [5] Raukema et al. used T s — 400 K [10j. 
For larger incidence angles they also reported data measured at a lower T s = 220 K. 
There, the lowering of the surface temperature led to an increase of the slope of the 
sticking curve above Ei = 1.0 eV [lOj . Unfortunately, no such lower temperature 
data was presented for normal incidence. If we assume that the trend seen at larger 
incidence angles prevails, we would expect normal incidence measurements at lower 
T s to yield a larger slope than the one of the experimental T s = 400 K curve shown in 
figure [5] - in closer agreement to the calculated sticking data. 

In this respect it is also important to note that from the observed non-monotonous 
dependence of So(E{) on T s and on the incidence conditions (Ei, incidence angle) 
Raukema et al. proposed the existence of two distinct dissociation mechanisms [TO] : 

(i) a direct one that governs dissociation at normal incident energies Ei > 1 eV and 

(ii) a precursor-mediated one at work in the threshold energy regime (Ei < 1 eV) , in 
which a molecular chemisorption state acts as the precursor. 

In the absence of substrate mobility and concomitant phononic dissipation channels 
in the model, the closest we can get to a precursor-mediated mechanism within the 
present theoretical approach are the molecularly trapped events we indeed observe 
for Ei > 1.1 eV. The probability for such trapping is rather energy independent and 
about 10~ 6 — 10~ 5 . Even if we thus assume each of these trajectories to eventually 
lead to dissociation, at this magnitude this contribution only significantly affects the 
total sticking coefficient at the direct onset around Ei — 1.1 eV as shown in figure [5] 
The resulting kinked sticking curve then has a shape that is very reminiscent of 
the experimental curve measured by Raukema et al. [10] for an off-normal incidence 
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(3.1) 







f , Eo (E) = [Nift^oT 1 E exp -ft (y/E-y/E^ 



(3.2) 



Non-adiabatic effects during the dissociative adsorption of O2 at Ag(lll)? 15 

angle of 30° and low surface temperature T s = 220 K. Unfortunately, again for normal 
incidence the reported data does not extend to similarly low incidence energies and 
surface temperatures to resolve this kink. Vice versa, we can at present only speculate 
if an account of substrate mobility in the calculations would for T s = 400 K modify 
the present frozen surface total sticking curve in figure [5] (dissociation plus trapping) 
to such an extent to reach quantitative agreement with the shown experimental curve. 
Particularly because of our scepticism that the here computed molecularly bound 
state on the adiabatic PES corresponds to the experimentally reported molecular 
chemisorption state, we do not believe that this will be the case. From the dependence 
on surface temperature Raukema et al. estimate a trapping probability of their 
precursor molecular chemisorption state of the order of 10~ 4 at incidence energies 
as low as Ei — 0.5 eV [10] . At such low energies our computed molecular state can 
not be accessed at all. 

Our current interpretation is therefore instead that the discrepancy in the initial 
sticking coefficient below Ei = 1 eV is primarily caused by another molecular precursor 
state that is absent in the theoretical calculations and that is not a property of pristine 
Ag(lll). This state would then also correspond to the strongly bound state seen in the 
TDS experiments 48, 51, 52, 52] . This interpretation that a molecular precursor state 
associated to surface imperfections is predominantly responsible for the dissociation 
probabilities in the intermediate energy range below Ei — 1.0 eV has already been 
advocated by Rocca et al. |14j to rationalize the discrepancy in their measured sticking 
data and the one of Raukema et al. precisely in this energy range. Our analysis fully 
supports this view and places an upper limit of 10~ 7 to the dissociative sticking at 
ideal Ag(lll) for incidence energies up to at least Ej = 0.9 eV. In this situation, new 
experiments are self-evidently required to finally settle this cause. Particularly for 
normal incidence they should be performed for a wide range of incidence energies (as 
in Raukema's work) and for as low surface temperatures as possible (as to not getting 
riddled by residual gas contamination). 

Finally, we come back to the original motivation and the possibility that the 
discrepancy between experiment and theory could also be due to non-adiabatic effects. 
However, such effects that either involve sudden charge transfer [IB] or spin flip 
hindrance [1] tend to release the system in an excited state and would thus in general 
further lower the dissociative sticking, which is already too low in its present adiabatic 
form. Due to the specific characteristics of the adiabatic PES not even such a lowering 
seems likely though. In contrast to e.g. the 02/Al(lll) system [HE], the adiabatic 
dissociative sticking at Ag(lll) is governed by very late energy barriers larger than 
1.1 eV and at Z-distances very close to the surface. Potential spin- or charge transfer- 
induced hindrances instead occur earlier in the entrance channel at Z-distances around 
2-3 A Q] 12] US]. Even if such hindrances would give rise to additional barriers, those 
molecules that are able to surmount them would thus still be confronted with the 
adiabatic barriers when they approach closer to the surface and into the PES region 
where the spin quenching or charge transfer has then already taken place. A noticeable 
spin-selection rule induced lowering of the dissociative sticking as compared to our 
present adiabatic curve would therefore only occur, if the additional restrictions in 
the entrance channel were larger than those due to the late adiabatic barriers. This is 
not possible for energetic reasons, because the present adiabatic MEP barrier is with 
1.1 eV already higher than the gas-phase O2 singlet-triplet gap of 0.98 eV [IS] ■ Similar 
arguments hold for possible non-adiabaticity effects due to charge transfer limitations, 
considering that the electron affinity of O2 is 0.44 eV [18]. 
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Figure 6. (Colour online) Comparison of the computed adiabatic initial 
dissociative sticking coefficient of O2 on Ag(lll) (black full circles, copied from 
figure^ and on Ag(100) (red open circles, taken from |17| ) for normal incidence 
conditions. The insets show 2D cuts of the minimum energy barrier configurations 
leading to dissociation on Ag(100) (upper inset) and on Ag(lll) (lower inset). 

As such we reach the same conclusions with respect to non-adiabatic effects as for 
O2 dissociation at Ag(100) [IT] . Also there, the dissociative sticking was found to be 
governed by large and late barriers above 1 eV on the adiabatic PES. Correspondingly, 
the computed sticking coefficients at Ag(100) and at Ag(lll) are strikingly similar as 
shown in figure [6] The higher values at Ag(100) for lower incidence energies are due 
to the higher accessible configurational space as illustrated in figure [6] by the elbow 
plots along the MEP. Intriguingly, at Ag(100) experiments do confirm the absence of 
dissociative sticking at 0.9 eV [S], which (considering the analogy of the theoretical 
findings at the two Ag surfaces) indirectly gives further support to our assessment 
that surface imperfections cause the larger dissociative sticking around of 10~ 4 at this 
energy in the experiments by Raukema et al. on Ag(lll). 

4. Summary and conclusions 

In summary, we have employed a DFT-based divide and conquer approach to study 
the reactivity of O2 molecules impinging onto the Ag(lll) surface. The determined 
adiabatic dissociative adsorption probability is extremely low and only exceeds 10~ 7 
for incidence energies above 1.1 eV. Our analysis directly relates this low reactivity to 
large energy barriers in excess of 1.1 eV located very close to the surface. This excludes 
that electronically non-adiabatic effects, either due to spin-selection rules [IJ[2] or due 
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to charge transfer limitations [TSJ [TB] , contribute to the low sticking coefficient. The 
late adiabatic barriers exceed both the gas-phase O2 singlet-triplet gap and the O2 
electron affinity. Even if non-adiabatic effects would give rise to additional barriers 
in the entrance channel, those molecules that are able to surmount them would thus 
still be confronted with higher barriers when they approach the adiabatic PES region 
closer to the surface. As such, the sticking coefficient is insensitive to a possible non- 
adiabaticity in the 02-Ag(lll) interaction and other quantities will have to be found 
as useful signatures for such effects. 

The highly activated nature of O2 dissociation at Ag(lll) is in good agreement 
with the existing experimental data. Still, quantitatively, large deviations of more 
than two orders of magnitude are obtained at intermediate incidence energies around 
0.8-0.9 eV, when taking the particular data set from Raukema et al. [TU] as reference. 
Instead, our data is well consistent with the estimate of 9 x 10~ 7 for this energy 
range by Rocca et al. [H] , who had already suggested surface imperfections as reason 
for the higher dissociation probability in Raukema's data. The analysis of possible 
uncertainties in the present approach fully supports this point of view. At best, the 
high surface temperature employed in the experiments by Raukema et al. could be 
an alternative reason, and in this respect it would be intriguing to repeat the present 
adiabatic calculations with a model that includes some account of substrate mobility. 
However, more important to finally settle the cause would be new experiments at lower 
surface temperatures and specifically investigating the role of surface imperfections 
(like the step density) for the dissociative sticking of O2 at Ag(lll). 
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